Momentum distribution, vibrational dynamics and the potential of mean force in ice 
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By analyzing the momentum distribution obtained from path integral and phonon calculations we 
find that the protons in hexagonal ice experience an anisotropic quasi-harmonic effective potential 
with three distinct principal frequencies that reflect molecular orientation. Due to the importance of 
anisotropy, anharmonic features of the environment cannot be extracted from existing experimental 
distributions that involve the spherical average. The full directional distribution is required, and we 
give a theoretical prediction for this quantity that could be verified in future experiments. Within 
the quasi-harmonic context, anharmonicity in the ground state dynamics of the proton is substantial 
and has quantal origin, a finding that impacts the interpretation of several spectroscopies. 
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Investigating the impact of hydrogen (H) bonding on 
molecular properties is the focus of intense research, but 
even behavior as fundamental as the equilibrium dynam- 
ics of the protons participating in H bonds remains poorly 
understood. Proton dynamics is reflected in the momen- 
tum distribution probed by deep inelastic neutron scat- 
tering (DINS) [J. Recent DINS studies of H bonded sys- 
tems have made striking observations, such as the pres- 
ence of a secondary feature in the tail of the spherically 
averaged distribution in confined water [5] , and estimates 
of a surprisingly large quantum kinetic energy of the pro- 
ton in undercooled water [3J. The secondary feature was 
attributed to quantum tunneling between the two wells 
of an anharmonic ID potential [2]. It is not clear, how- 
ever, to what extent the dynamics of an interacting many 
body system can be reduced to that of a single proton 
along a bond. For instance, it has been pointed out 
that anisotropy can mimic features of a spherical dis- 
tribution that one might associate to anharmonicity in a 
ID model [I], and yet so far there is no conclusive study 
of this issue. To interpret experiments in confined and 
undercooled water, the unknown details of the molecular 
structure are a severe source of difficulty. However, even 
in the simpler case of ice Ih, it is not clear if the physics 
can be captured by simple model potentials, and how an- 
harmonicity, anisotropy and structural disorder influence 
the momentum distribution. 

In order to tackle these issues we consider the 
open path integral Car-Parrinello molecular dynamics 
(PICPMD) data for ice Ih that yielded the accurate 
spherical momentum distribution reported in a prior pub- 
lication 5 . In this prior study, no attempt was made 
to relate the distribution to the equilibrium dynamics of 
the proton or to investigate the role of the environment 
in terms of a potential of mean force. In simulations this 
task is facilitated by access to the full 3D distribution, 
in contrast to experiments on polycrystalline samples, 



where only the spherically averaged distribution could 
be measured [U [6]. In addition, crystalline symmetry 
allows the use of harmonic analysis to quantify the rela- 
tion between the momentum distribution and vibrational 
dynamics, thereby elucidating the role of anharmonicity 
and disorder on the proton ground state. 

We find that anisotropy stemming from the molecu- 
lar orientations in the crystal has a larger effect on the 
momentum distribution than anharmonicity. The latter 
is effectively described within a quasi-harmonic model 
and is particularly important in the stretching motions, 
corroborating pump-probe laser experiments on the ex- 
cited state dynamics of ice and water [71 IH]- This finding 
impacts the interpretation of infrared and x-ray spectro- 
scopies, and regarding DINS experiments, the large effect 
of molecular anisotropy implies that it is not possible to 
unambiguously attribute to anharmonicity features of the 
spherically averaged distribution. Substantially more in- 
formation, capable of disentangling anisotropy from an- 
harmonicity, can be extracted from the directional dis- 
tribution, for which we now present the first theoretical 
prediction for a realistic system. 

The simulations of Ref. 5 sampled the end-to-end dis- 
tribution of the open Feynman paths of the protons [9] , 
i.e. v{x) — ;^rX]i^( x ) where the sum runs over the 
N p protons in the cell and the vector x points from 
one end of the path to the other. The momentum dis- 
tribution v(p) is the Fourier transform of z?(x). For 
each distribution £j(x) we compute the correlation ma- 
trix Ci t0l p — (x a xp). Within the statistical errors of 
the simulation the eigenvalues {c 2 }^ =1 of C{ are the 
same for all the protons, while the associated eigenvec- 
tors {vi,k}k—i are proton specific directions related by 
crystalline symmetry to the directions of the other pro- 
tons. This suggests an anisotropic Gaussian form for the 
end-to-end distribution: t'j(x) oc exp (— ^x T C^ x). A 
quantile analysis reported in the supplement [10] fully 
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supports this hypothesis. Thus the momentum distri- 
bution is fi(p) oc exp (— 2p-p T Cip), implying that the 
corresponding potential of mean force has the effective 
harmonic form V(r) = k^r T 'A^r, where M and r denote 
the proton mass and position. Ai has eigenvalues urf, and 
shares with Cj the eigenvectors, v^. The are related 
to the <j\ by, 
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and Wfe and v,-^ are denoted the principal frequencies and 
directions of proton i. Since the principal frequencies do 
not depend on i all the protons have equivalent local 
environments within the simulation error bars. 

By averaging over the protons we obtain the frequen- 
cies u>k with error bars in the first row of Tabic [IJ In terms 
of the a\ the spherically averaged end-to-end distribution 
takes the form, 
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Fig. [ija) shows that this curve differs negligibly from 
the corresponding "raw" distribution extracted from the 
simulation, indicating that an effective harmonic model 
faithfully represents the spherically averaged data. Con- 
sistent with chemical intuition, the associated principal 
directions reflect the orientation of each water molecule 
in the crystal. The principal axes corresponding to the 
highest frequency are close to the oxygen-oxygen near- 
est neighbor directions, whereas the eigenvectors associ- 
ated with the middle and lowest frequency correspond re- 
spectively to directions in and perpendicular to the HOH 
molecular plane. 

The PICPMD principal frequencies differ from their 
harmonic counterparts (see Table [I]). The latter were 
obtained with the phonon calculation discussed below. 
Thus the model that better represents the data is 
anisotropic and quasi-harmonic. We can now resolve, 
in the case of ice, a major issue that troubled the inter- 
pretation of experiments [I] by quantifying the relative 
importance of anisotropy and anharmonicity. We depict 
in Fig. [I] (b) the spherical distributions corresponding 
to, respectively, the quasi-harmonic model (first row of 
Table [I]), the harmonic model (second row of Table 0, 
and the isotropic model with frequency Q = 1186 cm -1 
that best fits the data. Anisotropy and anharmonic- 
ity are both significant, but anisotropy clearly has the 
larger effect. The isotropic model corresponds to a clas- 
sical Maxwell-Boltzmann distribution with an effective 
temperature T = 869i"T. In spite of T being signifi- 
cantly higher than the equilibrium temperature of the 
simulation (T — 269K), the isotropic model severely un- 
derestimates quantum effects, a finding that is also il- 
lustrated by a kinetic energy (Ek — lllmeV) approx- 
imately 30 percent smaller than the simulation value 
(E K = 143meV). 





wi(cm ) 


w 2 (cm 1 ) 


tj 3 (cm _1 ) 


E K {meV) 


PICPMD 


2639 ± 60 


1164 ± 25 


775 ± 20 


143 ±2 


Harmonic 


3017.6 ±8.2 


1172.5 ±8.9 


870.3 ± 14.6 


157.5 ±0.3 



TABLE I: Average proton principal frequencies and kinetic 
energies obtained from PICPMD and phonon calculations. 
The error bars reflect statistical errors and physical effect of 
disorder in the PICMD and phonon data, respectively. 
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FIG. 1: (color online) (a) The spherical end-to-end distribu- 
tion directly collected from PICPMD data (red dashed line) 
compared with that reconstructed by the anisotropic fit (blue 
line), (b) Comparison of the spherical momentum distri- 
bution of the harmonic crystal (black dot-dashed line) with 
anisotropic (blue line) and isotropic (red dashed line) fits. 



All the principal frequencies in Table [T] are well in ex- 
cess of the equilibrium temperature, indicating that the 
dynamics of the proton is dominated by quantum zero- 
point motion. Dependence of the molecular orientations 
upon the crystalline framework originates anisotropics 
that reflect the symmetry of the environment in the mo- 
mentum and end-to-end distributions. To study these 
effects we focus on the latter distribution, which factor- 
izes into the product of a spherical free-particle contri- 
bution and an anisotropic environmental component ny , 

^ Mk B Tx 2 ^ 

i.e. n(x) cx e ^f? ny(x) |llj . Rather than extracting 
ny(x) directly from the PICPMD data, which would be 
affected by substantial noise, we reconstruct ?V( X ) from 
the superposition of the individual proton contributions 
within the quasi-harmonic model. Here we use the fact 
that there are 24 unique orientations of the molecules in 
the hexagonal ice crystal [IS], and we also include the 
effects of proton disorder estimated below in the phonon 
calculation. Fig. [2] (a) depicts the log scale plot of one in- 
dividual environmental end-to-end distribution projected 
on the basal plane of ice Ih. The elliptic shape of the 
contour comes directly from the quasi-harmonic model. 
Fig. [2] (b) illustrates the log scale plot of the superpo- 
sition of all the environmental end-to-end distributions. 
The hexagonal shape of superpositioned distribution is a 
striking manifestation of quantum mechanics as in classi- 
cal physics ny (x) is equal to 1 . While the distribution is 
spherical at the center, hexagonal character emerges at 
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intermediate displacements and becomes pronounced in 
the tail of the distribution where blurring of the contour 
lines due to disorder can be detected. Experiments on 
ice Ih have only measured the spherical distribution [5J 
but it is likely that the full three dimensional distribution 
should become accessible in the future with improved in- 
strumentation and preparation techniques. Directional 
momentum distributions have already been reported for 
materials such as KDP [T3j and Rb 3 H(S0 4 ) 2 [H]. It 
should be noted, however, that the greatest sensitivity to 
anisotropy is in the exponential tail of the distribution, 
a finding indicating that substantial resolution may be 
necessary to experimentally disentangle anisotropy, an- 
harmonicity and other environmental effects. 




FIG. 2: (color online) (a) "Environmental part" of the end-to- 
end distribution corresponding to one individual proton pro- 
jected in the basal plane of ice Ih plotted in logarithmic scale, 
(b) "Environmental part" of the end-to-end distribution cor- 
responding to the superposition of all protons projected in 
the basal plane of ice Ih plotted in logarithmic scale. The su- 
per positioned end-to-end distribution reflects the symmetry 
of the oxygen sub-lattice. The blurring of the contour lines 
reflects the disorder effect detected in the phonon calculation. 

Lastly, we discuss the relationship between the princi- 
pal frequencies and the vibrational spectrum. The latter 
includes four main features experimentally: a stretch- 
ing band centered at « 3250 cm -1 [15], a bending 
band centered at « 1650 cm -1 [15], a wide librational 
band between ps 400cm -1 and 1050cm -1 [T5][T7] and a 
band of network modes below s» 400cm" 1 [T5]. These 
features are reproduced in the phonon spectrum of ice 
that we calculate by diagonalizing the dynamical ma- 
trix. This calculation is performed with Qbox [TU] by 
adopting the same supercell, electronic structure parame- 
ters and disordered proton configuration of the PICPMD 
simulation [5]. The dynamical matrix is calculated with 
a finite difference method (grid size of 0.0053A). The 
resulting phonon density of states shown in Fig. [3] (a), 
agrees with experiment, and is consistent with previous 
calculations [25], which did not include proton disorder, 
indicating that such effects have a small influence on the 
spectrum. We indicate phonon frequencies and eigenvec- 
tors by L)fr. h and ei ayk , respectively, where a are Carte- 
sian components, i, k = 1, • • • ,3N — 3, and N is the 
number of supercell atoms. In the quantum harmonic 
approximation the momentum distribution of particle i 



of mass Mi has the anisotropic Gaussian form fi(pi) oc 
exp ( — t;pJ Cf h Pi) with correlation matrix [2Tj . 
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FIG. 3: (color online) (a) Density of states of the phonon 
spectrum, (b) The population function for the principal axes 
corresponding to u>i (blue dot-dashed line), u)2 (red solid line) 
and u>3 (black dashed line). Network modes below 500cm -1 
contribute non-negligibly to all principal frequencies. 



As a consequence of disorder the eigenvalues of Cf a m , 
depend on the proton index i. The harmonic average 
frequencies are reported in the second row of Table [I] 
The corresponding standard deviations originate almost 
entirely from ice disorder, being at least an order of mag- 
nitude larger than the numerical errors estimated from 
the small asymmetry of the calculated dynamical ma- 
trix. The statistical errors in the PICPMD simulation 
(Table [T]) are on average a few times larger than the har- 
monic estimate of disorder, confirming that, within error 
bars, all proton environments are equivalent. We expect 
that longer runs combined with better estimators of the 
end-to-end distribution [IT] should improve the statis- 
tical accuracy to the point that disorder effects could 
become measurable in future simulations. 

The population function, 
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gives the weight of the phonon k in the principal direction 
I and is depicted in Fig. [3](b). It is found that Q\ is 94% 
stretching, u>2 is 47% bending and 48% libration, and W3 
is 97% libration. Taking only stretching, bending, and 
libration into account, and using weights proportional to 
h we infer that uii ~ 3160cm -1 , H12 ~ 1210cm -1 , and 
Q 3 r~j 895cm -1 . In comparison, the values in the second 



4 



line of Table [I] are red-shifted by contributions from net- 
work modes (6%, 4%, and 3% to U!\,u>2, an d 0J3, respec- 
tively), an intriguing result suggesting that fine details 
of the momentum distribution should reflect intermedi- 
ate range order properties of the H bond network. 

The potential energy surface is generated with the 
same protocol in path integral and phonon calculations. 
We thus attribute the difference between the average 
principal frequencies in the two rows of Table [T] to an- 
harmonicity. This originates from quantum derealiza- 
tion, present in the PICPMD simulation, which causes 
the proton to sample the potential over an extended 
range. Along the bond direction the proton spans from 

— 0.2A to ss +0.3A relative to the centroid of the path. 
This is much larger than the corresponding classical ther- 
mal spread (« ±0.05A) indicating that quantum anhar- 
monicity is essentially unaffected by temperature. The 
asymmetry of the quantal spread suggests that the first 
correction to the harmonic potential depends cubically 
on displacement. A visualization of the approximate ef- 
fective potential along the bond based on further calcula- 
tions is given in the supplement [10] and shows that while 
cubic terms dominate in the ground state, higher or- 
der corrections become important at displacements larger 
than ss 0.3A. In the supplement we also report harmonic 
estimates of the quantum effects on oxygen, which are 
non-negligible albeit smaller than for the protons. 

In conclusion, we find that to a large extent the mo- 
mentum distribution in ice is a simple anisotropic Gaus- 
sian distribution. This does not mean, however, that 
ice behaves like a harmonic crystal as the principal fre- 
quencies of the distribution differ from those of a har- 
monic crystal. Anharmonicity, enhanced by H bonding, 
is appreciable in the libration dominated Q3 and is par- 
ticularly significant in the stretching dominated u)\, in 
agreement with optical pump-probe experiments [71 [8] . 
The quantal character of the anharmonicity is consistent 
with the observed T-independence of the lifetime of ex- 
cited stretching modes in ice [7J. Our findings have im- 
plications for the calculation of observables in ice, such 
as infrared spectra, which typically ignore quantum an- 
harmonicity [22j . and x-ray absorption spectra, which 
typically ignore quantum configurational disorder |23j . 
The approach presented here could be applied directly 
to the study of other crystalline H bonded systems, and 
is also an important step towards a better understand- 
ing of the proton momentum distribution in disordered H 
bonded systems such as water under different conditions. 
In such cases only the spherically averaged momentum 
distribution is accessible in experiment and simulation 
can provide essential microscopic information to supple- 
ment and interpret the experimental data. Finally, we 
remark that while the qualitative picture emerging from 
our calculations is robust, the path integral data have 
relatively large error bars and the quantitative details 
depend on the accuracy of the underlying Born Oppen- 



heimer potential energy surface. The latter should reflect 
the known limitations of the GGA functional used in this 
study [24l [25] and comparisons with future high resolu- 
tion experiments should help to clarify this issue. 
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